{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "54652539",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from numpy import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "3babed15",
   "metadata": {},
   "outputs": [],
   "source": [
    "class HMM:\n",
    "    def __init__(self,A,B,pi):\n",
    "        self.A =A\n",
    "        self.B =B\n",
    "        self.pi =pi\n",
    "    def _forward(self,obs_seq):\n",
    "        # 取A = N x N\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        F = np.zeros((N,T))\n",
    "\n",
    "        # alpha = pi*b\n",
    "        F[:,0] = self.pi *self.B[:,obs_seq[0]]\n",
    "\n",
    "        for t in range(1,T):\n",
    "            for n in range(N):\n",
    "                # 计算第t时，第n个状态的前向概率\n",
    "                F[n,t] = np.dot(F[:,t-1], (self.A[:,n])) * self.B[n, obs_seq[t]]\n",
    "        return F\n",
    "    def _backward(self,obs_seq):\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        X = np.zeros((N,T))\n",
    "        # 表示X矩阵的最后一列\n",
    "        X[:,-1:] =1\n",
    "\n",
    "        for t in reversed(range(T-1)):\n",
    "            for n in range(N):\n",
    "                # 边权值为a_ji\n",
    "                X[n,t] = np.sum(X[:,t+1] * self.A[n,:] * self.B[:,obs_seq[t+1]])\n",
    "\n",
    "        return X\n",
    "    def viterbi(self,obs_seq):\n",
    "\n",
    "        N = self.A.shape[0]\n",
    "        T = len(obs_seq)\n",
    "\n",
    "        prev = np.zeros((T-1,N),dtype = int)\n",
    "\n",
    "        V = np.zeros((N,T))\n",
    "        V[:,0] = self.pi * self.B[:,obs_seq[0]]\n",
    "\n",
    "        for  t in range(1,T):\n",
    "            for n in range(N):\n",
    "                # 计算delta(j)*a_ji\n",
    "                seq_probs = V[:,t-1] * self.A[:,n] * self.B[n,obs_seq[t]]\n",
    "                # 记录最大状态转移过程\n",
    "                prev[t-1,n] = np.argmax(seq_probs)\n",
    "                V[n,t] = np.max(seq_probs)\n",
    "        return V,prev\n",
    "\n",
    "    def build_viterbi_path(self,prev,last_state):\n",
    "        \"\"\"\n",
    "        returns a state path ending in last_state in reverse order.\n",
    "        \"\"\"\n",
    "        T = len(prev)\n",
    "        yield(last_state)\n",
    "        # 从T-1开始，每次下降1\n",
    "        for i in range(T-1,-1,-1):\n",
    "            yield(prev[i,last_state])\n",
    "            last_state = prev[i,last_state]\n",
    "\n",
    "    def state_path(self,obs_seq):\n",
    "        V,prev = self.viterbi(obs_seq)\n",
    "\n",
    "        # build state path with greatest probability\n",
    "        last_state = np.argmax(V[:,-1])\n",
    "        path = list(self.build_viterbi_path(prev,last_state))\n",
    "\n",
    "        return V[last_state,-1],reversed(path)\n",
    "    \n",
    "    def baum_welch_train(self, observations, criterion=0.05):\n",
    "        n_states = self.A.shape[0]\n",
    "        # 观察序列的长度T\n",
    "        n_samples = len(observations)\n",
    "\n",
    "        done = False\n",
    "        while not done:\n",
    "            # alpha_t(i) = P(o_1,o_2,...,o_t,q_t = s_i | hmm)\n",
    "            # Initialize alpha\n",
    "            # 获得所有前向传播节点值 alpha_t(i)\n",
    "            alpha = self._forward(observations)\n",
    "\n",
    "            # beta_t(i) = P(o_t+1,o_t+2,...,o_T | q_t = s_i , hmm)\n",
    "            # Initialize beta\n",
    "            # 获得所有后向传播节点值 beta_t(i)\n",
    "            beta = self._backward(observations)\n",
    "\n",
    "            # 计算 xi_t(i,j) -> xi(i,j,t)\n",
    "            xi = np.zeros((n_states, n_states, n_samples - 1))\n",
    "            # 在每个时刻\n",
    "            for t in range(n_samples - 1):\n",
    "                # 计算P(O | hmm)\n",
    "                denom = np.sum(alpha[:, -1])\n",
    "                for i in range(n_states):\n",
    "                    # numer[1,:] = 行向量，alpha[i,t]=实数，slef.A[i,:] = 行向量\n",
    "                    # self.B[:,observations[t+1]].T = 行向量,beta[:,t+1].T = 行向量\n",
    "                    numer = alpha[i, t] * self.A[i, :] * self.B[:, observations[t + 1]].T * beta[:, t + 1].T\n",
    "                    xi[i, :, t] = numer / denom\n",
    "\n",
    "                # 计算gamma_t(i) 就是对j进行求和\n",
    "                gamma = np.sum(xi, axis=1)\n",
    "                # need final gamma elements for new B\n",
    "                prod = (alpha[:, n_samples - 1] * beta[:, n_samples - 1]).reshape((-1, 1))\n",
    "                # 合并T时刻的节点\n",
    "                gamma = np.hstack((gamma, prod / sum(prod)))\n",
    "                # 列向量\n",
    "                newpi = gamma[:, 0]\n",
    "                newA = np.sum(xi, 2) / np.sum(gamma[:, :-1], axis=1).reshape((-1, 1))\n",
    "                newB = copy(self.B)\n",
    "\n",
    "                # 观测状态数\n",
    "                num_levels = self.B.shape[1]\n",
    "                sumgamma = np.sum(gamma, axis=1)\n",
    "                for lev in range(num_levels):\n",
    "                    mask = observations == lev\n",
    "                    newB[:, lev] = np.sum(gamma[:, mask], axis=1) / sumgamma\n",
    "\n",
    "                if np.max(np.abs(self.pi - newpi)) < criterion and \\\n",
    "                                np.max(np.abs(self.A - newA)) < criterion and \\\n",
    "                                np.max(np.abs(self.B - newB)) < criterion:\n",
    "                    done = True\n",
    "                self.A[:], self.B[:], self.pi[:] = newA, newB, newpi\n",
    "    def simulate(self, T):\n",
    "        def draw_from(probs):\n",
    "            # np.random.multinomial 为多项式分布，1为实验次数，类似于投掷一枚骰子，丢出去是几，probs每个点数的概率，均为1/6\n",
    "            # 给定行向量的概率，投掷次数为1次，寻找投掷的点数\n",
    "            return where(np.random.multinomial(1, probs) == 1)[0][0]\n",
    "\n",
    "        observations = np.zeros(T, dtype=int)\n",
    "        states = np.zeros(T, dtype=int)\n",
    "        states[0] = draw_from(self.pi)\n",
    "        observations[0] = draw_from(self.B[states[0], :])\n",
    "        for t in range(1, T):\n",
    "            states[t] = draw_from(self.A[states[t - 1], :])\n",
    "            observations[t] = draw_from(self.B[states[t], :])\n",
    "        return observations, states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "id": "3ac27aab",
   "metadata": {},
   "outputs": [],
   "source": [
    "states =('sunny','cloudy','rainy')\n",
    "\n",
    "observations = ('dry','dryish','damp','soggy')\n",
    "\n",
    "start_probability ={'sunny':0.6,'cloudy':0.2,'rainy':0.2}\n",
    "\n",
    "transition_probability = {\n",
    "    'sunny':{'sunny':0.5,'cloudy':0.375,'rainy':0.125},\n",
    "    'cloudy':{'sunny':0.25,'cloudy':0.125,'rainy':0.625},\n",
    "    'rainy':{'sunny':0.25,'cloudy':0.375,'rainy':0.375}\n",
    "}\n",
    "\n",
    "emission_probability ={\n",
    "    'sunny':{'dry':0.6,'dryish':0.2,'damp':0.15,'soggy':0.05},\n",
    "    'cloudy':{'dry':0.25,'dryish':0.25,'damp':0.25,'soggy':0.25},\n",
    "    'rainy':{'dry':0.05,'dryish':0.1,'damp':0.35,'soggy':0.50}\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "id": "562440d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5   0.375 0.125]\n",
      " [0.25  0.125 0.625]\n",
      " [0.25  0.375 0.375]]\n",
      "[[0.6  0.2  0.15 0.05]\n",
      " [0.25 0.25 0.25 0.25]\n",
      " [0.05 0.1  0.35 0.5 ]]\n",
      "[0, 1, 2, 3]\n",
      "[0.6 0.2 0.2]\n",
      "forward: P(O|lambda) = 0.026441\n",
      "backward: P(O|lambda) = 0.026441\n"
     ]
    }
   ],
   "source": [
    "def gengerate_index_map(labels):\n",
    "    index_label = {}\n",
    "    label_index = {}\n",
    "\n",
    "    i =0\n",
    "    for l in labels:\n",
    "        index_label[i] =l\n",
    "        label_index[l] =i\n",
    "\n",
    "        i+=1\n",
    "    return label_index,index_label\n",
    "\n",
    "states_label_index,states_index_label = gengerate_index_map(states)\n",
    "observations_label_index,observations_index_label = gengerate_index_map(observations)\n",
    "\n",
    "def convert_observations_to_index(observations,label_index):\n",
    "    list =[]\n",
    "    for o in observations:\n",
    "        list.append(label_index[o])\n",
    "    return list\n",
    "\n",
    "\n",
    "def convert_map_to_matrix(map,label_index1,label_index2):\n",
    "    m = empty((len(label_index1),len(label_index2)),dtype = float)\n",
    "    for line in map:\n",
    "        for col in map[line]:\n",
    "            m[label_index1[line]][label_index2[col]] = map[line][col]\n",
    "    return m\n",
    "\n",
    "def convert_map_to_vector(map,label_index):\n",
    "    v = empty(len(map),dtype = float)\n",
    "    for e in map:\n",
    "        v[label_index[e]] =map[e]\n",
    "    return v\n",
    "\n",
    "A = convert_map_to_matrix(transition_probability,states_label_index,states_label_index)\n",
    "print (A)\n",
    "B = convert_map_to_matrix(emission_probability,states_label_index,observations_label_index)\n",
    "print (B)\n",
    "observations_index = convert_observations_to_index(observations,observations_label_index)\n",
    "print (observations_index)\n",
    "pi = convert_map_to_vector(start_probability,states_label_index)\n",
    "print (pi)\n",
    "\n",
    "h = HMM(A,B,pi)\n",
    "# 人为定义的海藻状态序列\n",
    "obs_seq = ('dry','damp','soggy')\n",
    "obs_seq_index = convert_observations_to_index(obs_seq,observations_label_index)\n",
    "\n",
    "# 计算P(o|lambda)\n",
    "F = h._forward(obs_seq_index)\n",
    "print (\"forward: P(O|lambda) = %f\" %sum(F[:,-1]))\n",
    "X = h._backward(obs_seq_index)\n",
    "print (\"backward: P(O|lambda) = %f\" %sum(X[:,0]*pi*B[:,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "a1d2373e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "最有可能的隐藏序列为： ['sunny', 'cloudy', 'rainy']\n",
      "viterbi: P(I|O) =0.010547\n"
     ]
    }
   ],
   "source": [
    "# 计算P(I|o)\n",
    "p,ss = h.state_path(obs_seq_index)\n",
    "path = []\n",
    "for s in ss:\n",
    "    path.append(states_index_label[s])\n",
    "\n",
    "print(\"最有可能的隐藏序列为：\" ,path)\n",
    "print(\"viterbi: P(I|O) =%f\"% p)\n",
    "\n",
    "# result：\n",
    "# 输入观察序列为：['dry','damp','soggy']\n",
    "# 最有可能的隐藏序列为： ['sunny', 'cloudy', 'rainy']\n",
    "# viterbi: P(I|O) =0.010547"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "id": "b5dc1b61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.28\n"
     ]
    }
   ],
   "source": [
    "# 参数估计\n",
    "observations_data, states_data = h.simulate(100)\n",
    "guess = HMM(array([[0.33, 0.33, 0.34],\n",
    "                       [0.33, 0.33, 0.34],\n",
    "                       [0.33, 0.33, 0.34]]),\n",
    "                array([[0.25, 0.25, 0.25, 0.25],\n",
    "                       [0.25, 0.25, 0.25, 0.25],\n",
    "                       [0.25, 0.25, 0.25, 0.25]]),\n",
    "                array([0.7, 0.15, 0.15])\n",
    "                )\n",
    "guess.baum_welch_train(observations_data)\n",
    "# 预测问题\n",
    "states_out = guess.state_path(observations_data)[1]\n",
    "p = 0.0\n",
    "for s in states_data:\n",
    "    if next(states_out) == s: p += 1\n",
    "\n",
    "print(p / len(states_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "fee67987",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 1 3 0 1 1 0 2 1 2 2 1 0 0 0 2 2 3 2 0 2 1 1 2 2 1 1 3 2 0 3 2 3 2 0 1 0\n",
      " 2 3 3 1 3 0 2 3 1 3 1 0 3 0 3 0 1 3 2 3 3 2 3 1 0 3 3 3 3 0 0 2 2 2 2 1 0\n",
      " 2 2 3 3 0 0 0 0 2 0 0 0 0 2 1 0 3 0 2 0 2 3 0 0 2 0]\n"
     ]
    }
   ],
   "source": [
    "print(observations_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "304994f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "source": [
    "print(type(observations_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "id": "68098a2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "a=np.array(\n",
    "    [\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/4+1/10.,1/4.-1/10.,1/4.,1/4.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "c4209aa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#定义观察序列\n",
    "q={'研': 0, '究': 1, '生': 2, '物': 3,'很': 4,'有': 5,'意': 6,\n",
    "   '思': 7, '他': 8, '大': 9, '学':10,'时':11,'代':12,'是':13,\n",
    "   '的':14, '专':15, '业':16, '首':17,'选':18,'目':19,'标':20}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "id": "cb816fb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def observation_list(ob):\n",
    "    ol=[]\n",
    "    T=len(ob)\n",
    "    for i in range(T):\n",
    "        ol.append(q[ob[i]])\n",
    "    return np.array(ol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "a71531b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3 4 5 6 7]\n"
     ]
    }
   ],
   "source": [
    "print(observation_list('研究生物很有意思'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "f3f51950",
   "metadata": {},
   "outputs": [],
   "source": [
    "test01=HMM(a,b,pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "id": "5ce85898",
   "metadata": {},
   "outputs": [],
   "source": [
    "test01.baum_welch_train(observation_list('研究生物很有意思'),criterion=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "id": "0914378b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.20938628 0.29963899 0.24548736 0.24548736]\n",
      " [0.20651574 0.30314743 0.24516842 0.24516842]\n",
      " [0.208      0.30133333 0.24533333 0.24533333]\n",
      " [0.208      0.30133333 0.24533333 0.24533333]]\n",
      "[[0.10671748 0.         0.         0.         0.         0.\n",
      "  0.81705574 0.07622677 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.04870624 0.         0.         0.         0.         0.\n",
      "  0.87011669 0.08117707 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.07862408 0.         0.         0.         0.         0.\n",
      "  0.84275184 0.07862408 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.07862408 0.         0.         0.         0.         0.\n",
      "  0.84275184 0.07862408 0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]]\n",
      "[0.35 0.15 0.25 0.25]\n"
     ]
    }
   ],
   "source": [
    "print(test01.A)\n",
    "print(test01.B)\n",
    "print(test01.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0772c72",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
